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ABSTRACT 

We present new observations of the CB130 region, composed of three separate 
cores. Using the Spitzer Space Telescope we detected a Class and a Class II 
object in one of these, CB 130-1. The observed photometric data from Spitzer and 
ground-based telescopes are used to establish the physical parameters of the Class 
object. SED fitting with a radiative transfer model shows that the luminosity 
of the Class object is 0.14 — 0.16 L©, which is a low luminosity for a protostellar 
object. In order to constrain the chemical characteristics of the core having the 
low luminosity object, we compare our molecular line observations to models 
of lines including abundance variations. We tested both ad hoc step function 
abundance models and a series of self-consistent chemical evolution models. In 
the chemical evolution models, we consider a continuous accretion model and an 
episodic accretion model to explore how variable luminosity affects the chemistry. 
The step function abundance models can match observed lines reasonably well. 
The best fitting chemical evolution model requires episodic accretion and the 
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formation of CO2 ice from CO ice during the low luminosity periods. This process 
removes C from the gas phase, providing a much improved fit to the observed 
gas-phase molecular lines and the CO2 ice absorption feature. Based on the 
chemical model result, the low luminosity of CB130-1 is explained better as a 
quiescent stage between episodic accretion bursts rather than being at the first 
hydrostatic core stage. 



Introduction 



The Spitzer Space Telescope Legacy Project, From Molecular Cores to Planet Forming 
Disks (c2d, lEvans et al. 2003) has completed a survey of nearby star-forming regions. It 



covered five clouds, containing 1024 objects classified as young stellar Objects (YSOs) (Evans 



et al. 2009 ). When a dense core forms a central protostar, it develops a luminosity source from 
the accretion of infalling mass. The accretion luminosity is given as Lace = GM^Macc/R*, 
where is the mass of the protostar. Mace is the mass accretion rate onto the protostar, 
and i?* is the radius of the protostar. According to the standard model, the mass accretion 
rate from spherical infall of the envelope is Maec — 2 x 10^® Mqjt^'^ if the infall occurs at 
the thermal sound speed at 10 K (Shu 1977). With a typical protostellar radius of 3 R©, 
and a mass at the stellar/brown dwarf boundary of 0.08 Mq, the resulting Laec is 1.6 Lq. 
Any luminosity from contraction of the forming star will add to this accretion luminosity, 
making it a minimum value in the standard model. 



In contrast to the predictions, the c2d survey showed that 59% of the 112 embedded 
protostars (Class O/I) have luminosity lower than 1.6 Lq ( Evans et al.||2009" ), indicating that 
either is very small, or Mace is lower than expected, or both. This luminosity problem 



was aggravated by the discovery of Very Low Luminosity Objects (VeLLOs, di Francesco 



etaL][2007| . VeLLOs are defined as embedded protostars with internal luminosity lower than 
0.1 Lq. The internal luminosity of an embedded protostar is the luminosity of the protostar 
and disk, excluding luminosity from external heating. By using the correlation between 70 



fim flux and the internal luminosity of protostars, Dunham et al. (2008) identify 15 VeLLO 



candidates in the full c2d sample. Several VeLLOs including L1014-IRS ( Young et al.|2004"a ), 
L1148-IRS (iKauffmann et al.||2005|, L1521F-IRS (iBourke et all [20061), IRAM 04191-IRS 



(Andre et al. 1999 Dunham et al. 2006), L328-IRS (Lee et al. 2009), and L673-7 (Dunham 



et al. 2010a), have been studied in detail. Among those VeLLOs, IRAM 0419H-1522 and 



L673-7 drive strong molecular outflows (Andre et al. 1999 Dunham et al. 2006, 2010a). A 



low accretion luminosity with a significant molecular outflow can imply higher Mace in the 



past (Dunham et al. 2010a). 
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One possible explanation for the sources with low luminosities combined with strong 
molecular outflows is that the mass accretion is not a constant process but rather episodic, 



and the VeLLOs are in quiescent phases between the mass accretion bursts (Kenyon et al. 



1990[ [Enoch et al.|2009[ [Dunham et al.|2010b[ ) . [Dunham etal] ( |2010b| ) present a set of evolu- 
tionary models describing collapse including episodic accretion. The luminosity distribution 
of YSOs can be matched only when the model includes episodic accretion. 

While dust continuum emission can provide the physical structure and evolutionary 
stage of a core, molecular line observations trace the dynamics and chemistry of the core. 
The simple empirical models of step function ( Lee et aT][2003 ) and drop function (J0rgensen 



et al. 2004 ) abundance proflles have been used to approximate the abundance profiles. The 



chemical evolution model during protostellar collapse developed by Lee et al. (2004) has been 



tested for individual cores, like B335 ( [Evans et al.|[2005| , and L43 ( [Chen et al.||2009| ). Using 



this model, Lee (2007) studied chemical evolution in VeLLOs. The low luminosity sources 



can provide chemical laboratories to test the astrochemistry of low luminosity environments. 
Because the chemical equilibrium time can be longer than the duration of an accretion burst, 
the abundance profile may provide a fossil record of previous luminosity bursts. 

This paper presents new observations of the CB130 region. With a detailed analysis of 
Spitzer, submm continuum, and molecular line data we will reveal the physical properties of 
this source and use CB130 as a laboratory to study the chemistry in low luminosity sources. 
In ^ we give a general introduction to CB130. The observations are described in ^ ^ 
presents images, photometry, and molecular line data, ^presents dust continuum radiation 
transfer models used to determine physical parameters. The line modeling with an empirical 
step function is described in ^ The chemical evolution models with different luminosity 
evolution are presented in ^ and we summarize our findings in ^ 



2. CB130 



CB130-1 first appears in the catalog of Clemens & Barvainis (1988), which presents 248 



small isolated molecular clouds. Clemens & Barvainis (1988) called the core CB130. Later 



Lee & Myers (1999) identified 3 cores near CB130 and named them CB130-1, CB130-2, and 



CB130-3, respectively, from south to north. Fig. [T] shows the R-band image from the Digital 
Sky Survey (DSS), with CS (left) and N2H'*" (right) molecular line contours (De Vries et 
al., in prep, see ® overlaid. The three cores are clearly identified as separate dark cores in 



the region. Wu et al. (2007) presented Submillimeter High Angular Resolution Camera II 



(SHARC-II) 350 iJ,m and 450 fim detections of CB130-1. [Dunham et al] ( [20081 ) identified an 
embedded source in CB130-1 with Lj„f ~ 0.07 Lq based on a correlation between 70 fim 
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flux and the internal source luminosity. However, this result is an estimate, and is only good 
to within a factor of ~ 2. Recently, Launhardt et al. (2010) presented a survey of low mass 



star forming cores in 32 Bok globules, including CB130-1. They found a faint near infrared 
(NIR) source and a very red star in CB 130-1, and they identified these as Class and Class 
II objects, respectively. 

CB 130 is located in the Aquila Rift cloud, whose galactic coordinates are 20° < / < 40° 
and —6° < b < +14°. The local standard of rest velocity of the Aquila Rift cloud is 



8 km/s (Dame & Thaddeus 1985), which is similar to the velocity of CB130-1 (see §6 



Straizys et al. (2003) determined the distance to the Aquila Rift, using extinction methods. 



to extend between 225 and 310 pc, and the central region at 270 pc. We thus adopt a 
distance to CB130-1 of 270 pc but note the substantial uncertainty. 



Observations 



We obtained observational data from various telescopes. All three instruments mounted 



on Spitzer observed CB130-1 in the c2d program. The Infrared Array Camera (IRAC; Fazio 



et al.||20O4 ), the Infrared Spectrograph (IRS; Houck et al.||2004 ), and the Multiband Imaging 



Photometer for SIRTF (MIPS; Rieke et al. 2004). The IRAC images were obtained on 



September 2nd 2004 and September 3rd 2004 (Program ID [PID] 139; c2d [Evans et~aL][2003 



AOR keys 0005145856 and 0005146368), the IRS spectra were obtained on April 26th 2006 
(PID 179; PI: N. J. Evans; AOR key 0015921153), and the MIPS 24 and 70 /im data were 
obtained on September 25th 2004 (PID 139, AOR key 0009412608, 0009421824). 

The IRAC and MIPS images from c2d were reduced by the Spitzer Science Center (SSC) 
using the standard pipeline version S13 to produce Basic Calibrated Data (BCD) images. 
The BCD images are obtained by subtracting the dark and bias levels, and performing flat 
fielding and sky subtraction. A full explanation of this extra processing data can be found 



in the documentation of the final delivery of data from the c2d legacy project (Evans et al. 



2007). 



The IRS data were reduced following the IRS pipeline version S14.0.0, which produces 
BCD files. IRS BCD files are two dimensional spectra, which have been processed including 
saturation flagging, dark-current subtraction, linearity correction, cosmic ray correction, 
ramp integration, droop correction, stray light removal or crosstalk correction, and fiat- 
field correction. Detailed explanations about the data delivery can be found in the c2d 



Spectroscopy Explanatory Supplement (Lahuis et al. 2006) 



Spitzer IRAC and MIPS 24 /im images of CB130-1 were also obtained in the Spitzer 
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GO-2 program (cores2deeper; PID 20386; PI: P. C. Myers), which obtained deeper images of 
a sample of cores from the c2d program. IRAC 4 band images were obtained on September 
16th 2005 (AOR key 0014606848) and on September 19th 2005 (AOR key 0014607104). 
MIPS 24 /im image was obtained on October 5th 2005 (AOR key 0014615040). Additional 
Spitzer MIPS images at 24 /im and 70 fim were obtained along with 160 yum image on 
May 18th 2007 (PID 30384; PI: T. L. Bourke; AOR key 0018159616). IRAC images were 
reduced in the same way as the c2d data. MIPS images from PID 20386 and PID 30384 



were reduced using the MIPS Data Analysis Tool (DAT; Gordon et al. (2005)). The data 
reduction process is described in Stutz et ah] ( |2007[ ). 



Submm Common User Bolometer Array (SCUBA)[^ observations at 850 yum, taken at 
the James Clerk Maxwell Telescope (JCMT), were obtained from the CADC archive. The 
SCUBA jiggle maps were analyzed using the standard SCUBA User Reduction Facility 
(SURF) reduction package and aperture photometry was calibrated using Uranus observa- 
tions obtained from the same night (see Shirley et al.|2000 for the flux calibration procedure). 
The beam size at 850 is 15". 

Near-infrared observations of CB130-1 were obtained during 24-25 March 2005 using 
the Infrared Sideport Imager (ISPI) on the 4-meter Blanco telescope at Cerro Tololo Inter- 
american Observatory (CTIO). Deep observations at J and H were obtained with 60-second 
and 30-second exposures, respectively. Deep Ks observations were obtained with 2 coadded 
10-second exposures. The total exposure times achieved for these deep observations were 9, 
76, and 17 minutes at J, H, and Ks, respectively. To obtain photometry of most moderately 
bright stars, which were saturated in these deep observations, shallower 3-second exposures 
were obtained. The individual images were reduced using IRAF and IDL procedures fol- 
lowing the standard method, as described in Huard et al. (2006), to construct the final 
stacked images at J, H, and Ks. Sources were identified in these final images and l'.'2-radius 



aperture photometry derived using Phot Vis, version 1.09 (Gutermuth et al. 2004), an IDL 



GUI-based aperture photometry tool that uses standard DAOPHOT packages (Landsman 



1993). Finally, the photometry was calibrated by comparison with published photometry 
of sources listed in the Two Micron All Sky Survey (2MASS) catalog. The average seeing 
during the J observations was ~ 1'.'3, and ~ I'.'l during the H and Ks observations. 



The molecular line spectral data were obtained from several telescopes (see Table |3] 
Isj). The majority of the lines were observed at the Caltech Submillimeter Observatory 



-'^This research used the facihties of the Canadian Astronomy Data Centre operated by the National 
Research Council of Canada with the support of the Canadian Space Agency. 
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(CSO) 0from 2005 to 2009. For the data obtained between 2005 - 2007, we used an acousto- 
optic spectrometer (AOS), with 1024 channels and 50 MHz bandwidth. For the 2008 - 
2009 data, we used a fast fourier transform spectrometer (FFTS) with 8192 channels and 
500 MHz bandwidth. We used the 230 GHz heterodyne receiver for all the lines except 
H2D"'" {Jk^iK+i = li,o for which we used the 345 GHz receiver. The position 

switched scans were taken when the optical depth at 225 GHz was between 0.1 and 0.2. 
The calibration uncertainty is about 10%. N2H"'"(J = 1 — )■ 0) and CS (J = 2 — > 1) maps 
were obtained at the Five College Radio Astronomy Observatory (FCRAO; De Vries et al., 
in prep) and N2H+ (J = 1 — )■ 0) data were obtained from the Arizona Radio Observatory 
(ARO) on March 2007. We reduced the molecular line data from the CSO, FCRAO, and 
ARO 12 m telescope using CLASS Q 



4. Results 

4.1. The Images and Photometry 

A three-color image of CB130-1 using cores2deeper IRAC data is presented in Fig. [2} 
with 3.6 fim as blue, 4.5 /xm as green, and 8.0 fim as red. Two YSOs are identified in these 



data based on their positions in various color-color & color- magnitude diagrams (Harvey et al. 



2007). These two YSOs are 15" apart, which corresponds to 4100 AU at a distance of 270 pc 
dStraizys et al.||2003|). We name the western source (a = 18^ 16"^ 16.4^ 6 = -02° 32' 38'.'0 



[J2000]) CB130-1-IRS1 and the eastern source {a = 18^* 16"^ 17.4^ 6 = -02° 32' 41'.'1 [J2000]) 
CB130-1-IRS2 and label these sources in Fig. |2] 

Fig. |3] shows a three color image of CB 130-1 using the NIR data with J as blue, H 
as green and Ks as red. Both YSOs are detected and labeled in the figure. Cone shaped 
nebulosity extends from CB130-1-IRS1 toward CB130-1-IRS2, as marked by arrows in the 
figure. Similar nebulosity is also observed in L1014 (Huard et al. 2006[ ) and is indirect 



evidence of an outflow. However, we mapped the CB130-1 region with CO (J = 2 — )■ 1) at 
the CSO and found no significant evidence of outflowing gas. 

Fig. |4] shows 850 fim contours overlaid on the 1.25 — 24 fim NIR and Spitzer images. 
The 850 /im data probes the envelope since dust emission from the envelope becomes optically 
thin at long wavelengths. The 850 /im contours peak at CB130-1-IRS1 and avoid CB130-1- 
IRS2. CB130-I-IRS2 is brighter at all wavelengths up to 8 /xm, but CB130-1-IRS1 is stronger 



^The Caltech Submillimeter Observatory is supported by the NSF. 
3http://www.iram.fr/IRAMFR/GILDAS 
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at 24 yum. These facts support the conclusion by Launhardt et aL (2010) that CB130-1-IRS1 
is younger and more embedded than CB130-1-IRS2. 

A color-color diagram constructed from the IRAC photometry for all sources detected 
in the Spitzer images of the CB130 region is plotted in Fig. [5j Sources inside the box are 
generally Class II sources ( [Alien et al.||2004| ). The two sources in CB130-1 are clearly redder 
than background stars. The color-color diagram shows that CB130-1-IRS1 and CB130-1- 
IRS2 are in different phases of evolution. CB130-1-IRS2 is clearly in the region of Class 
II sources. The sources for which either the [3.6]- [4.5] color is greater than 0.8 mag or the 
[5. 8] -[8.0] color is greater than 1.1 mag are all faint and classified as galaxies except for 



CB130-1-IRS1 (Harvey et al. 2007). Since no YSO candidates in CB130-2 and CB130-3 



were detected with CTIO and Spitzer, we will concentrate on CB130-1 for the remainder of 
the paper. 



4.2. CB130-1-IRS1 

The photometric data for CB130-1-IRS1 are given in Table [T| which lists wavelength, 
flux density, uncertainty in flux density, aperture diameter, and telescope. At 3.6 — 24 
fim, we list the flux averaged between the c2d and cores2deeper images, weighted by the 
inverse of the uncertainties. The SED is plotted in the upper panel of Fig. [6} including the 
IRS spectrum, which clearly contains 10 fim silicate and 15.2 fim CO2 absorption features. 
Using the 15.2 /im CO2 absorption feature, we can derive the CO2 abundance toward CB130- 
1-IRSl. The continuum for the spectrum is constructed by fitting a first order polynomial 
to the ranges 12 — 14.7 and 16.2 — 18 /im. After subtracting the continuum, we calculated 
the CO2 column density from A^(C02)= / r(z/)(ii//y4, where r(t/) is an optical depth, u is 



the wave number (cm ^), and A = 1.1 x 10 is the band strength (Gerakines et al. 1995). 



We find / T{u)du = 28 ± 0.4 cm'^ and N{C02)= (2.6 ± 0.05) x 10^^ cm'^. We calculated 



the H2 column density using 850 /im data. Using equation (6) in Lee et al. (2003), we have 
iV(H2)= 4.5 X 10^2 ^^-2 p^Qj^ A^(C02) and A^(H2), the abundance of CO2 ice is 5.8 x IQ-^ 
which requires a substantial fraction of the carbon. 

Based on the observed SED, we can calculate standard observational signatures in- 
cluding Lboi, Thoi, and i^feo«/-^(>35o/im)- We calculate Lj,oi by integrating the flux over all 
wavelengths using the average of the results obtained from midpoint and prismoidal integra- 
tion methods. is the temperature of a blackbody having the same fiux-weighted mean 



frequency as the source and is calculated following Myers & Ladd (1993). ^(>35o/im) is the 
luminosity emitted at A > 350 /xm. CB130-1-IRS1 has 1^ = 0.23 ± 0.03 Lq, = 59 ± 1.6 
K, and -^vfeo;/L(>35o/im) = 9.1. Lboi/L(^^3^o^m) is used to distinguish Class from Class I 
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sources. The pre-protostellar core/Class transition occurs around -f'6oz/-^(>35o/im) ~ 35 
and the Class O/I transition occurs at -Z^6o«/-^(>35o/im) ~ 175 (Young & Evans 2005). Since 
Lboi/ L (^^350 fim) is less than 35, we would guess that CB130-1 is in the pre-protostellar core 
stage from these criteria, but CB130 clearly contains embedded sources. Another criterion 
to define a Class object is the bolometric temperature of the source. Class and I objects 
are divided at T^oi = 70 K (Chen et al. 1995); thus CB130-1-IRS1 is classified as a Class 
source based on Tboi- 

The bolometric luminosity includes the luminosities of the embedded object, the disk 
and the ISRF, which can add several tenths of a Lq to Lboi (Evans et al. 2001). The 
luminosity of the internal source can be estimated in two ways. First, assuming that all 
emission at A > 100 fim arises from external heating, L^nt can be calculated by integrating 
the fiux at A < 100 /im, giving Lint ~ 0.056 Lq. Alternatively, we can estimate Lint by 
applying the relation between the fiux at 70 yum, normalized to 140 pc, and Lint found by 
Dunham et al. (2008), giving Lint ~ 0.07 L0. These two estimates indicate that CB130-1- 



IRSl is in any case a low luminosity object and may be a VeLLO. However, they are only 
rough estimates; radiative transfer models of CB 130-1 are required to reliably determine the 
physical properties of the internal source (^. 



Assuming that the core is optically thin at submillimeter wavelengths, the envelope 
mass can be estimated from M = -]^^^^-, where 5*,^ is the fiux density at 850 yum, D is the 
distance to the object, By{Td) is the black-body function, and is the opacity of gas and 



dust at 850 /xm (Young et al. 2003). Pontoppidan's dust model (Pontoppidan et al., 2010, 



in preparation) gives Ky 



1.02 X 10-2 cm-i 



-1 



and Ky = 1.8 x 10 



cm 



dust (Ossenkopf & Henning 1994). If we assume the core is isothermal at 



-1 for 0H5 
10 K, the 



estimated envelope mass within a 120" diameter aperture (radius of 16,200 AU) is 3.5 Mq 
with the Pontoppidan model, and 2.0 Mq with the 0H5 model. These estimates will be 
refined later. 



4.3. CB130-1-IRS2 



The SED of CB130-1-IRS2 is plotted in the lower panel of Fig. [6j CB130-1-IRS2 is not 
detected at A > 24/im. The slope of the SED of CB130-1-IRS2 from 2.17 /im to 24 /xm is 



—0.7, which indicates that CB130-1-IRS2 is a Class II object (Greene et al. 1994). We can 



estimate the bolometric luminosity of the source by integrating the SED. The bolometric 
luminosity of CB130-1-IRS2 is Lboi = 0.071 ± 0.003 L0 and the bolometric temperature is 
Tboi = 1150 ± 15 K. These values of Lboi and Tboi are not corrected for extinction along the 
line of sight, and are thus lower limits to the intrinsic values. 
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4.4. The Molecular Line Data 

Molecular line observations can trace the dynamics in the core, such as infall, rotation, 
and outflow. The observed molecular lines for CB130-1, CB130-2, CB130-3 are summarized 
in Table [3} Table [4], and Table |5] respectively. The detected molecular lines are plotted in 
Fig. [7] as solid lines. Toward CB130-2, the detected molecular lines are plotted in the upper 
three panels in Fig. [8j Toward CB 130-3, the detected lines are plotted in the lower three 
panels in Fig. |8| 

In the next section, we will focus on CB130-1-IRS1, since the source is a deeply em- 
bedded Class object. We will flnd the internal luminosity and chemical properties of the 
source. 



5. Radiative Transfer Model 

As discussed above, a detailed radiative transfer model is required to accurately de- 
termine the internal luminosity of CB130-1-IRS1. The internal luminosity will determine 
whether the source is a VeLLO or not. This section describes one-dimensional and two- 
dimensional models of this source. These models trace the radiation from a source in a 
dusty region, including scattering, absorption, and re-emission by the dust and generate 
model SEDs that can be compared to the observations. 

The one-dimensional models provide one estimate of the internal luminosity and a spher- 
ical model that is needed for analyzing the molecular lines. The two-dimensional models 
provide an alternative estimate of the internal luminosity and a better flt to the SED, at the 
expense of having more free parameters. 



5.1. One-dimensional radiative transfer models 

We modeled CB130-1-IRS1 with the one-dimensional radiative transfer code, DUSTY 



(Ivezic et al. 1999). DUSTY calculates the dust temperature proflle and SED of a protostar 
embedded in a dense core. The parameters of the model divide into those associated with 
the envelope, or core, and those associated with the internal source. 

The envelope is characterized by the outer radius, inner radius, density distribution, 
dust properties, and external radiation fleld. The outer radius of the envelope is not well 
constrained by the radiative transfer modeling as long as it is large enough to include all 
the emission. We flx the outer radius at 35000 AU, based on the angular size of CB130-1 
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(3'0 X 1^3, Lee & Myers 1999). The inner radius [Rin) will be a free parameter. We assume 
a power-law radial density distribution for the envelope (n = nQ{r /tq)'"), where tiq is the 
density at a fiducial radius, tq. The value of no was constrained to match the long- wavelength 
emission, and a = 1.8 ±0.2 was constrained to match the radial intensity profile. See Shirley 



et al. (2002) for an explanation of this method. 



We use the Pontoppidan dust opacity model (Pontoppidan et al., 2010, in preparation). 



The envelope is heated by both the embedded source and the ISRF described by Evans 



et al. (2001). The ISRF can vary in strength and is modified by extinction by a surrounding 



low density envelope. To constrain the ISRF, we used an initial model of the envelope and 
modeled the CO (J = 2 — > 1) line, since the CO (J = 2 — > 1) line is sensitive to the UV 
part of the ISRF, and thus the extinction. The best fit value for the UV part of ISRF is 
Go = 4.5 X 10"'^, corresponding to Ay = 3. We will discuss this in ^ With these dust 
properties, the density normalization (no) yields the total mass of the model envelope within 
35,000 AU to be 5.3 Mq. This mass depends only slightly on the best-fitting parameters for 
luminosity and inner radius, as long as they are close to the best-fit values. Within 16,200 



AU (the size used to estimate the observed mass in section 4.2), the model contains a mass 
of 2.2 M0, slightly less than the isothermal calculation within that radius. 

The internal source is characterized by an internal luminosity (Lint) and a stellar tem- 
perature [Ts). The luminosity of the disk is included in the model, but the disk luminosity 
is not easily separable from the stellar luminosity, so we treat the sum as the free parameter: 
Lint = Ls + Ldisk- The emission from the disk is averaged over all inclination angles and 



added to the input spectrum (Butner et al. 1994). The 24 /xm photometry data constrains 



the disk parameters. The surface density profile of the model disk follows S(r) oc with 
p = 1.5. We choose a disk temperature distribution following T(r) oc r^^. When the disk 
is a fiat passive disk, q = 0.75. For the fiat passive disk, all the luminosity comes from the 
re-radiation from the central star. A fiat passive disk produces less 24 /xm emission than 
observed (lower panel of Fig. [9]) . We varied the index of the temperature distribution until 
it fitted the 24 /xm observed data. The best fit power law index of the disk temperature 
distribution is g = 0.4, which corresponds to a fiared disk (Chiang & Goldreich 1997). The 



disk temperature distribution is also less steep than the fiat disk case if the disk has an in- 
trinsic luminosity source, such as accretion and back warming from the envelope in addition 
to radiation from the central star. 

DUSTY calculates the radiative transfer from the internal source and the external radi- 



ation field. Then, the OBSSPHere (Shirley et al. 2002) program, which convolves the model 



with a beam, is applied to the DUSTY result to convert models into observed fiux densities 
and model radial profiles at wavelength longer than 100 /xm. The model parameters are 
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constrained by the observed SED and the observed radial intensity profile at 850 /im. 

We tested a grid of models to obtain the best-fit free parameters. After the constraints 
described above, the remaining free parameters are Tg, Lint, and We calculated the 
reduced to quantify the quality of the fit. The data used to constrain the model covers 
3.5 < A < 850 /im, with a total of nine points of photometry data. We do not include 
H band photometry data, because H band flux comes mainly from scattered photons in the 
envelope, and DUSTY modeling does not properly treat scattering at short wavelengths. 
Also, we focus on the photometric data, and check the IRS spectrum only for consistency. 
We have 9 data points and 3 free parameters, so the degree of freedom of reduced is 6. 



The X plots are in Fig. 10 When we plot Fig. 10 , we hold fixed the other two parameters 
at the values with the smallest x^ values. We plot in panel (a) of Fig. 10 varying from 
lOOOK to 7000 K. The value of does not critically affect the modeled SED, since all 
the stellar photons are reprocessed by dust close to the star. The lowest is found at 
Ts = 3000 ± 500 K. Panel (b) is obtained by varying L^nt = + L^isk- The best fit 
is found at Lint = 0.16 ± 0.005 L©. The L^nt is a little higher than the VeLLO criterion, 
but much lower than the accretion luminosity calculated from the mass accretion rate in 
Shu (1977)'s model. Panel (c) is the x^ plot varying The inner radius of the envelope 



controls the optical depth. When the optical depth increases, the flux density decreases at 
short wavelengths because dust absorbs energy at short wavelengths and re-emits at long 
wavelengths. The best fit inner radius is = 350 ± 25 AU, which gives r = 0.22 at 100 
//m. 

The resulting model SED is plotted in Fig. |9| The flux density at wavelengths longer 
than 100 /xm comes mainly from the envelope. Since the envelope is bigger than the aperture 
size, the modeled flux density is higher than the observations. To compare to observations, we 
used OBSSPHere to convolve the model with the aperture that was used for the photometry. 
The convolved model data are marked with open circles in Fig. [9j The SED shows a good fit 
to the observed data at 24 yum, 70 /xm, 160 /xm, 350 yum, and 850 yum. However, the modeled 
SED shows a 24 — 30 /im turnover and shallow 15.2 /im CO2 feature, which is not an ideal 
fit to the IRS spectrum. 



5.2. Two-dimensional radiative transfer models 

The best fit result with the one-dimensional radiative transfer model fails to fit the IRS 
spectrum at wavelength 24 — 30 /im. The flux density in the mid-infrared region mainly 
comes from the disk emission. There are limitations to models of disks averaged over all 
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inclinations. For a more realistic explanation of the MIR spectrum, we used the two di- 



mensional axisymmetric Monte Carlo code, RADMC (DuUemond & Dominik 2004). The 



RADMC program solves the radiative transfer equation in 2-D circumstellar dust configu- 
rations. Once a dust temperature distribution is determined, RAYTRACE is used to create 
smooth SEDs for various inclination angles. As with the 1-D models, we use some constraints 
to determine some parameters before running a fine grid; this reduces the multi-dimensional 
fitting problem to a manageable size. 

We used the same dust opacity and ISRF models used for the 1-D model. For the 
2-D model envelope, we used a slowly rotating cloud core with a small angular momentum 



(Terebey, Shu, & Cassen 1984, hereafter TSC). The initial state of the TSC model is a 
singular isothermal sphere, with density proportional to r~^. As the core evolves, the collapse 
of the isothermal sphere includes the effects of an initially uniform and slow rotation. The 
density of the TSC envelope is determined by the initial angular velocity, i^o and the time, t. 
We set f2o = 10~^^s~^, which is the value used in TSC model. As t increases, the infall radius 
gets larger and the density distribution gets flattened. The best fit density distribution is 
constrained, as described for the 1-D model, to correspond to t = 35000 yr. We use the inner 
radius and the outer radius as 350 AU and 35000 AU respectively, which were determined by 
the 1-D model and the core size. The mass of the envelope that matches the long-wavelength 
emission is 5.5 Mq, only slightly larger than the 1-D model estimate in §5.1. 



For the disk density structure, we use the axisymmetric flared disk (DuUemond et al. 



2001 Whitney et al.|2003 Pontoppidan et al.|2007b ). The central disk structure is as follows: 

S(i?) 



p{R,z 

where the surface density S is 



exp 



2H,{R)\ 



1) 



and the disk scale height is given as 



-'disk 



R 



R. 



disk 



Hp{R) H( 



disk 



R 



R. 



■disk 



R 



R, 



■disk 



(2) 



(3) 



We set the disk outer radius as 350 AU, which is the inner radius of the envelope. The 
disk inner radius is 0.03 AU, set to the location at which a dust temperature is 1500 K. We 



use most of the other disk model parameters in Whitney et al. (2003). The flaring index afi 
is set to 0.25, and the surface density power law index is set to p = 1. The disk mass is set 
as Mdisk=0-Ol Mq. 
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We ran a grid of models varying three free parameters to fit the observed SED. As 
for the 1-D models, we calculated for the photometric points to determine the best fit. 
The free parameters are the disk scale height, Hdisk/Rdisk at i?rfjsfc=100 AU, the internal 
luminosity. Lint, and the inclination angle. Again, we plot the versus the free parameter. 



setting the other two free parameters at particular values. Panel (a) of Fig. 11 presents x 



,2 



varying the disk scale height Hdisk/Rdisk in the range between 0.01 (Whitney et al. 2003) 



and 0.125 (Pontoppidan et al. 2007b[ ). The best fit occurs for Hdisk/Rdisk = 0.05 ± 0.005. 



Panel (b) in Fig. 11 is the x^ vs internal luminosity. Since the disk is now heated only from 
the reprocessing of the central star, the internal luminosity affects the SED mostly in the 
MIR. If the internal luminosity is low, it cannot heat the disk enough. If Lint is high, the 
modeled SED becomes stronger than the observed flux density everywhere. The best flt is at 
Lint = 0.14 ± 0.005 Lq, which is a little bit lower than found for the one- dimensional model 
(0.16 Lq). Panel (c) is the x^ plot function of the inclination angle. The inclination 
angle affects the MIR SED, since the inclination angle affects the disk emission. There is 
a broad minimum in x^ around an inclination angle of 70°, but the SED for the model at 
70° badly underestimates the observations at 24 /im, which are well fitted by a model with 
inclination angle of 50°, which however badly overestimates the observations at 70 /im. The 



SEDs with these two inclination angles are plotted in Fig. [12] Neither model matches all of 
the IRS spectrum either. In particular, the observed CO2 feature is deeper than produced 
by any of the models. 

The best fit parameters are summarized in Table. |6} The internal luminosity from the 
one-dimensional model is 0.16 L©, the internal luminosity from the two-dimensional model 
is 0.14 L0. Both luminosities are slightly higher than the VeLLO criterion {Lint = 0.1 L0). 
CB130-1-IRS1 is a low luminosity source, but not a VeLLO. 



6. Molecular Line Model 

Chemical modeling of the line observations is needed to determine the chemical prop- 
erties of the source. We found that the chemical properties can give further clues to the 
evolutionary state of the source. In this section, we use step-function models of the chemical 
abundances, which will provide a simple comparison to the full chemical models employed 
in the next section. 

Because the chemical models are not very sensitive to the 2-D structure of the source, 
we use 1-D models of the dust. For each chemical model, we first calculate the dust radiative 
transfer to get the dust temperature distribution. Next, we calculate the gas kinetic temper- 
ature (§ 6.1), and then use step function abundance models, which can give a reasonable fit 
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to the observed molecular lines (Lee et al. 2003 Evans et al. 2005) (§ 6.2). A Monte Carlo 



code ( Choi et al.||1995 ) calculates the molecular excitation and a virtual telescope simulation 
produces line profiles, which can be compared to the observed lines (§ 6.3). Comparing 



observations to line profiles predicted from empirical abundance profiles can constrain the 
abundance profiles. 



6.1. Gas Temperature 



We calculate the radial gas temperature distribution using a gas energetics code (Doty 



& Neufeld 1997; Young et al. 2004b), which includes energy transfer between gas and dust. 



heating by cosmic rays, photoelectric heating and molecular cooling. 

The collision rate between gas and dust depends on the cumulative grain cross-section. 
In this simulation, the grain cross section per baryon is 6.09 x 10~^^ cm^ ( Young et al.|2004b ), 
based on the 0H5 model. This is slightly different from the opacity model we have used for 
DUSTY modeling, but the line modeling does not critically depend on the opacity model. 
We take the cosmic ray ionization rate as 3.0 x 10~^^ s~^ (Roberts & Herbst 2002; van der 



Tak & van Dishoeck 2000) 



The photoelectric heating is determined by the strength of UV part of the ISRF and 
the electron number density. The electron number density is assumed to be 0.001 cm~^ 
dDoty et al.||2002D . The CO line shows stron ger emission when the envelope is heated by the 
unattenuated ISRF. So we vary the extinction of the ISRF by the surrounding material until 



it reproduces the observed CO line. See Evans et al. (2005) for a full explanation of this 



method. The relevant part of the ISRF for the gas heating is in the FUV, where the ISRF 
is characterized by Go- Assuming a plane parallel slab, the ISRF is attenuated following 
Go = G{ISM)e-^uv ^here G{ISM) is the unattenuated FUV field. The best match to the 
observed CO line requires Gq/G{ISM) = 4.5 x 10~^. The relation between optical depth in 
the UV band and the V band is tuv = 1-8 Ay. Thus we attenuate the ISRF at the outer 
boundary of the envelope by Ay = — 0.56/n(4.5 x 10~^) = 3. We used this value for the 
external attenuation of the ISRF in §5.1[ Of course, the FUV is further attenuated as it 
propagates into the cloud, and the model accounts for this. 

The cooling rate mainly depends on the abundance and the line width of CO. We set 



the CO fractional abundance relative to H2 to be 7.4 x 10 ^ ( Bacmann et al.|[2002 J0rgensen 
et al. 2002). We use a Doppler b = 0.8 km s~^, which fits best the CO molecular line. 



The calculated gas temperature is plotted in the middle panel of Fig. 13, along with 
the dust temperature. The dust temperature and gas temperature are almost the same in 



- 15 - 



the inner region due to the high gas-dust colhsion rate because of the high density in the 
inner region, shown in the upper panel. In the outer part of the cloud the gas temperature is 
lower than the dust temperature, because the gas temperature has decoupled from the dust 
temperature, and the photoelectric heating from the attenuated ISRF is too low to heat the 
gas. 



6.2. Abundance Profile 

We use step function abundance profiles. Molecules tend to freeze out onto the grain 
surfaces in cold and dense regions. The inner part of the cloud has higher density compared 
to the outer region, so molecules freeze out faster. If the dust temperature becomes high 
enough, the molecules will sublimate. The CO sublimation temperature is about 20 K. The 
highest temperature in the model cloud is less than 20 K in the innermost region of the 
envelope, 0.003 pc, since CB130-1-IRS1 is a low luminosity source. The dust temperature 
is too low to sublimate the frozen-out CO. In the outermost part, the density is too low for 
significant freeze-out. Thus, a step function is a reasonable profile for the cloud. The three 
parameters of the abundance profile are the undepleted abundance (Xq), depletion radius 
(r/)), and depletion fraction (/d)- The best-fitting parameters are summarized in Table. [7j 



The abundance profiles are plotted in the last panel in Fig. 13 



6.3. Monte Carlo simulation 



The line profiles are modeled with the Monte Carlo (MC) code (Choi et al. 1995) to 



calculate level populations of molecules. The MC code can handle arbitrary 1-D distributions 
of the systematic velocity, the density, the kinetic temperature, the microturbulence, and 
the abundance self-consistently. 

Infalling envelopes sometimes exhibit a stronger blue wing for optically thick lines with 



a substantial excitation temperature gradient (Evans 1999; Lee et al. 2007). Even though 



CB130-1 contains a protostar, there is no significant feature of infall in molecular lines. All 
lines show a single peak, so we do not include any systematic velocity in the models. 

Doppler b parameters are determined by the width of each line. We have used b = 0.8 
km s~^ for CO. Since the CO line is optically thick and dominated by the outer part of the 
envelope, the CO line is difficult to fit with the same Doppler b parameter as others. We 
used b = 0.2 km s~^ for all the other lines, based on the widths of relatively optically thin 
lines, such as C^^O and H^^CO. 
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After MC modeling, a virtual telescope (VT) program is used to integrate the emission 
along the line of sight and to match the velocity and spatial resolutions, as well as the beam 
efficiency. All lines are assumed to be centered at 7.77 km s~^, determined by averaging the 
LSR velocity of optically thin lines. 

By comparing modeled spectra to observations, we can estimate the abundance. The 
best fit abundances are given in Table [7j We set the depletion radius as 0.04 pc. We used the 



isotope ratio for the low-mass protostellar envelope described in J0rgensen et al. (2004). The 



ratio of O to ^^O is 540, and the ratio between C and ^^C is 70. CO destroys N2H+ molecules, 
so N2H"*" has a reverse abundance profile compared to the CO abundance. The N2H"'" line 
has a hyperfine structure. The MC code calculates line modeling for each hyperfine line, and 
then they are added together. 

The best fit results are shown in Fig. (t) The modeled C^^O, CS, and N2H+ ( J = 1 -> 0) 
lines show reasonable fits to the observed lines, while modeled CO and HCO^ lines show a 
poor fit with the observed lines. The strengths are in a reasonable range for both HCO"*" 
and H^'^CO"'", but the abundance that fits the observed H^'^CO^ line produces a self-reversed 
feature in the modeled HCO^ line. The modeled HCO+ has a low excitation temperature 
in the outer region, so it shows absorption features at the center of the line. However, the 
observed HCO^ line does not show a self-reversed structure. The modeled H2CO also shows 
a self-reversed structure. The prediction of self- reversed lines that are not observed is a 



common problem (Carr et al. 1995 Chen et al. 2009), and it suggests that macroturbulent 
and clumpy 3-D models would be more appropriate. Though modeled DCO"*", CS, and N2H"'" 
(J = 3 — )■ 2) line strengths are weaker than the observed lines, the difference is less than a 
factor of 2.5. 



7. Chemo-dynamical Models 

7.1. Chemical Models with Luminosity Evolution Scenarios 

CB130-1 is a low luminosity source. There is no significant evidence of a CO outfiow, 
but there is a refiection nebula suggestive of an outfiow cavity. CB130-1 could be in the 
first hydrostatic core stage or in a quiescent stage between episodic accretion bursts. With 
chemical evolution modeling, we may be able to distinguish these possibilities. Because 
it takes time to relax the chemistry in the core, evidence of a past high luminosity stage 
might be present in the abundances. We use the evolutionary chemo-dynamical model by 



Lee et al. (2004) to test this idea. This model calculates the chemical evolution of a core 



from the prestellar core to the embedded protostellar core stage. At each time step, the 
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density profile, the dust temperature, the gas temperature and abundances are calculated 
self-consistently. 

At the prestellar stage, the density distribution is described by a sequence of Bonnor- 
Ebert spheres. We adopt the time duration at the preprotostellar core stage in Table 4 in 



Chen et al. (2009), which has a longer total timescale at the prestellar stage. Once a singular 
isothermal sphere has been reached with n(r) oc r~^, the inside-out collapse is initiated, and 
n{r) approaches r~^'^ inside the infall radius. The infall radius propagates outward at the 
sound speed through the envelope. 

We use three kinds of luminosity evolution models in the accretion stage. First, in 



Model 1, we adopt the accretion luminosity from Young & Evans (2005) including the First 



Hydrostatic Core stage (FHSC; preoE|[l969| [B^p93l [M^unaga et al.||1998| ). The FHSC 
represents the core in an initial quasi-static contraction stage when it is still molecular. When 
the core temperature reaches 2000 K, the molecular hydrogen dissociates, so the core has a 
second collapse and the FHSC stage ends. Models indicate that the FHSC has a short, but 



poorly determined, lifetime, about 10 years, (Boss & Yorke 1995 Masunaga et al. 1998) 



a very low luminosity, and a high optical depth. While two candidates have recently been 
proposed ( |Chen et al.||2010| [Enoch et al.]|2010D , the FHSC has not yet been clearly detected. 
In Model 1, the FHSC lasts 20000 years and the collapse of the FHSC occurs in 100 years. 
When t = 20000 yr. Lace ~ 0.011 Lq, which is only 10% of the internal luminosity of CB130- 
1-IRSl. The luminosity increases rapidly after 20000 years. We adopt a time at the very 
end of the FHSC stage of Model 1, which has a luminosity much lower than that of CB130- 
1-IRSl. To match the luminosity of the object, we would have to be catching the source in 
an extremely short-lived phase while the FHSC is collapsing. 



The luminosity evolution of Model 2 is taken from Lee et al. (2004). In Model 2, the 



FHSC lasts 20000 years and the transition occurs in 32000 years, providing a better chance 
to see a source with the observed luminosity. However, this longer transition stage has no 
particular theoretical justification. Here, we choose 45000 years as our desired time since 
Lace is 0.16 Lq at that time. 

The last model of luminosity evolution. Model 3, uses the episodic accretion luminosity 



model from Dunham et al. (2010b). The episodic accretion is included in a simple, idealized 



way rather than a self-consistent model. The model assumes no accretion from disk to star 
most of the time, but continued infall from envelope to disk. A mass accretion burst occurs 
when the disk mass reaches 0.2 times the stellar mass. Then mass accretion from the disk 
to protostar increases to M = 1 x 10~^ MqJt~'^, which is about 100 times the average mass 
accretion rate. The accretion luminosity depends on M, so the accretion luminosity also 
evolves episodically. The envelope mass of CB130-1 is 5.3 Mq— 5.5 Mq. So among the 
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three masses in Dunham et al. (2010b), a 3 Mq envelope model is the closest. However, 
the luminosity evolution of the model with a 3 Mq envelope never falls below 1 Lq after 
the FHSC stage ends. So we take the luminosity evolution of the 1 Mq envelope model in 



Dunham et al. (2010b). In Model 3, we adopt 78000 years as our desired time step. At 78000 



yr, the internal luminosity is 0.14 Lq. The cloud has gone through 6 burst episodes, and 
3000 years have passed after the sixth burst. The luminosity evolution of the three models 



are plotted in Fig. 14 The early stage evolution is the same in all three models. After 20000 
years. Model 2 has the lowest luminosity among the three models. 

We calculate the dust temperature distribution with DUSTY at each time step. We use 
the same parameters described in § |5.1[ Then we calculate the gas temperature considering 



the parameters described in § 6.1 The chemical evolution is calculated for each of 512 gas 



parcels as it falls into the central region. Gas inside the infall radius carries the memory 
of the conditions from farther out. The chemical calculation includes interactions between 
gas and dust grains and gas-phase reactions, but does not include chemical reactions on 
grain mantles explicitly. We assume the surface binding energy of species onto bare Si02 
dust grains. We use an updated chemical network with a new binding energy of N2 and 



deuterated species, which is used in Chen et al. (2009) 



The abundance profile at each time step is calculated through the chemical evolution 
model. We use the same isotope ratio as in § |6.3[ The abundance profiles are plotted in 



Fig. 15 The abundance profiles of Model 1 and Model 2 do not have a CO evaporation 
region at the center. In both cases, the central source has a low luminosity, and it does not 
heat the envelope enough to make CO sublimate. Model 2 has higher abundances compared 
to Model 1, except for DCO+. In Model 3, the abundances of CO and C^^O show the CO 
sublimation region at the central region because of past periods of high luminosity. The 
modeled abundances of HCO+, H^'^CO^ are all low compared to the best fit step function 
abundances. 



Once all the physical parameters are determined, we use MC code (Choi et al. 1995) 



to calculate the level populations. Then we use VT code to produce line profiles. The 
molecular line profiles with Model 1 are plotted in Fig. 16, and with Model 2 in Fig. 17 
The modeled lines of HCO"'", H^^CO^, CS, and N2H"'" are weaker than the observed lines, 
and the DCO"^ line is stronger than the observed line. N2H"'" shows a better fit with Model 
2, but the J = 3 — )■ 2 line is still weak. The molecular lines with Model 3 are plotted in 



Fig. 18 In Model 3, the C^^O line is much stronger than observed lines. On the other hand, 
HCO^, H^^CO"*", and N2H^ are weaker than observed lines. The modeled N2H^ (3"2) line 
is very weak. All three models show strong C^^O and weak N2H"'" lines compared to the 
observed lines. The modeled C^^O is much stronger in Model 3, because the episodic bursts 
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have sublimated the CO in the inner region. The N2H"'" hne is weak in all three models. 

All the models have caveats in their luminosity evolution scenarios. Model 1 does not 
actually match the observed Lint ~ 0.14 Lq— 0.16 Lq. Model 2 assumes that the transition 
stage of FHSC lasts for 32000 years, which is unrealistic. It can give a source luminosity 
of 0.16 Lq, and it gives a slightly better fit to the C^^O and N2H"'" compared to Model 1. 
However, the modeled DCO"^ line is too strong, and the modeled CS line is too weak. In 
Model 3, the envelope mass is lower than our estimates for CB130-1. Also, the best fit TSC 
model density from the 2-D model is 35000 years, while the desired luminosity of Model 3 
is reached at t = 78000 years, which means that more envelope material from Model 3 has 
already fallen in. The problems with the luminosity evolution scenarios show how chemical 
modeling can provide additional constraints that are distinct from those from modeling the 
continuum. 



7.2. Including the Conversion of CO into CO2 Ice 



We have one more trick to try. As noted in ^4.2, a substantial amount of carbon is 



contained in CO2 ice. Once CO freezes out onto dust grain surfaces, a fraction of CO 



ice turns into CO2 ice (Pontoppidan et al. 2007a). CO2 ice does not sublimate during the 



accretion burst even when the CO evaporation temperature is reached. So the amount of CO 
returned to the gas phase is less than the amount of CO frozen onto grain surfaces during the 
quiescent phase. The net effect is to decrease gas phase CO, and enhance N2H"'". We used 
the same episodic accretion luminosity evolution model as Model 3, and varied the fraction 
of CO ice which turns into CO2 ice in the chemical network (Model 4). When 80% of CO ice 
turns into CO2 ice, the modeled C^^O and N2H+ lines match well with the observed lines. 
Also DCO"^ and CS fit the data better than in Model 3. The abundance of this model is 



plotted in Fig. 15 and the molecular line profiles are plotted in Fig. 19 The CO2 ice column 
density obtained from the best fit model is 2.94 x 10^^ cm~^, 1.1 times observed value, within 
the likely range of systematic uncertainties. 

Among the models, the episodic model including CO2 ice formation fits best to the 
observed lines in all molecules. Based on the modeling result, we suggest that CB130-1 
has gone through episodic accretion, forming CO2 ice out of CO ice in quiescent phases. 
This model can explain the low luminosity of the source, the strong CO2 ice feature, and 
most of the gas phase emission lines. There are still inconsistencies in this model, so further 
development of episodic accretion models is needed. 
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8. Summary 

We presented a detailed study of a low luminosity object in the CB130-1 region. We 
performed radiative transfer modeling and chemical modeling to explain the core with a low 
luminosity protostar in it. 

The embedded protostar CB130-1-IRS1 has Lint = 0.14 L0 (1-D model) to Lint = 0.16 
L0 (2-D model). This is a slightly higher luminosity than a VeLLO has, but still CB 130-1- 
IRSl is a low luminosity source at about 0.1 of the expected luminosity for steady accretion 
onto an object at the boundary between stars and brown dwarfs. 

We tested both a step function abundance model and self-consistent chemical evolution 
models. The step function abundance profile fits observed lines reasonably well, but is of 
course, ad hoc. For the self-consistent models, we use three different luminosity evolution 
scenarios to explain the low luminosity of the object. All three luminosity evolution models 
with the standard chemical network show that the modeled C^^O fine is strong and N2H+ 
is weak compared to the observations. The step function model has more free parameters 
than chemical evolution model, so it is more feasible to fit the observed lines. 

The deep 15.2 //m CO2 ice feature indicates that CO2 ice has formed from CO ice. We 
added a reaction to the chemical network, which turns CO ice into CO2 ice. A model with 
that reaction decreases the C^^O abundance and increases the N2H+ abundance. With the 
episodic accretion model and the modified network, we can explain the low luminosity of 
the source, the strong CO2 ice feature, and most of the gas phase emission lines at the same 
time. So CB 130-1 is most likely neither a FHSC nor a very low mass object, but a more 
evolved protostar in a quiescent stage between accretion bursts. With a further study of 
CO2, CO and N2H"'" in low luminosity objects, we may find chemical imprints of episodic 
accretion in low luminosity objects. 
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Table 1. Photometry of CB130-1-IRS1 



A 


5.(A) 




o 


Aperture diameter 


Telescope 


(/iin) ( 


mJy) ( 


mJy) 


(arcsec) 




1.25 < 0.736 








CTIO 


1.64 


0.088 


0.002 


2 


CTIO 


2.15 


< 1.27 








CTIO 


3.6 


1.58 




0.08 


1 7 




4.5 


4.57 




0.23 


1.7 


Spitzer 


5.8 


6.12 




0.31 


1.9 


Spitzer 


8.0 


7.28 




0.36 


2.0 


Spitzer 


24.0 


52.4 




5.2 


6.0 


Spitzer 


70.0 


312 




32 


50 


Spitzer 


160.0 


5840 


1700 


100 


Spitzer 


350.0 


2800 




700 


40 


cso 


850.0 


1480 




296 


120 


JCMT 


Note. - 


- The 350 


(Um continuum SHARC-II data is pre- 


sented by 


Wu et al. 


( 


2007 


)■ 




Table 2. 


Photometry of CB130-1-IRS2 


A 


5.(A) 


a 


Aperture diameter 


Telescope 


(/xm) 


l^mJy) (m 


Jy) 


(arcsec) 




1.25 


2.46 0.005 


2 


CTIO 


1.64 


7.59 0.007 


2 


CTIO 


2.15 


12.7 




0.0 


2 


CTIO 


3.6 


17.6 




0.9 


1.7 


Spitzer 


4.5 


18.2 




0.9 


1.7 


Spitzer 


5.8 


17.8 




0.9 


1.9 


Spitzer 


8.0 


20.5 




1.0 


2.0 


Spitzer 


24.0 


27.8 




2.8 


6.0 


Spitzer 
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Table 3. Observation parameters and molecular line data for CB130-1 



Line 


Frequency 


Vmb 


Beam size 






Av 


n 


Telescope 




(GHz) 




(arcsec) 


(K km s- 


(km s^-"- 


) (km s-l) 


(K) 




CO 2-1 


230.537970 


0.7 


33 


9.38 


8.27 


4.73 


1.86 


CSO 


Ci^'O 2-1 


224.714368 


0.8 


34 








< 0.14 


CSO 


Ci*0 2-1 


219.560352 


0.8 


35 


0.48 


7.61 


0.71 


0.64 


CSO 


HCO+ 3-2 


267.557620 


0.7 


29 


0.82 


7.67 


0.60 


1.3 


CSO 


DCO+ 3-2 


216.112605 


0.7 


36 


0.26 


7.68 


0.47 


0.53 


CSO 


H"CO+ 3-2 


260.255339 


0.8 


30 


0.09 


7.78 


0.49 


0.17 


CSO 


N2H+ 3-2 


279.511701 


0.7 


27 


0.45 


7.62 


0.88 


0.10 


CSO 


N2H+ 1-0 


93.173700 


0.6 


70 


2.27 


7.35 


1.07 


0.99 


ARO 


N2H+ 1-0 


93.176258 


0.5 


54 


1.30 


7.43 


0.44 


0.38 


FCRAO 


CS 2-1 


97.980953 


0.5 


54 


0.29 


7.78 


0.54 


0.43 


FCRAO 


N2D+ 3-2 


231.321635 


0.7 


33 








< 0.24 


CSO 


H2D+ 1,1,0 - 1,1,1 372.421340 


0.7 


23 








< 0.29 


CSO 


H2CO312 — 2ii 


225.697787 


0.7 


34 


0.21 


7.74 


0.53 


0.37 


CSO 



Table 4. Observing parameters and molecular line data for CB 130-2. 



Line 


Frequency rjmb Beam size 


/ TXdv vlsr 


Av 




Telescope 




(GHz) 


(arcsec) 


(K km s ^) (km 


s-i) (kl 


tn s"^) 


(K) 




CO 2-1 


230.537970 0.7 


33 


10.1 


7.99 


5.54 


1.72 


CSO 


C^^O 2-1 


219.560352 0.8 


35 


0.40 


7.19 


0.56 


1.04 


CSO 


HCO+ 3-2 267.557620 0.7 


29 


0.23 


7.23 


0.47 


0.47 


CSO 


N2H+ 3-2 


279.511701 0.7 


27 








< 0.13 


CSO 


N2D+ 3-2 


231.321635 0.7 


33 








< 0.30 


CSO 


N2H+ 1-0 


93.176258 0.5 


54 








< 0.24 


FCRAO 


CS 2-1 


97.980953 0.5 


54 


0.32 


7.38 


0.52 


0.58 


FCRAO 
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Table 5. Observing parameters and molecular line data for CB130-3. 



Line 


Frequency rj^b Beam size 


jTXdv 


vlsr 






Telescope 




(GHz) 


(arcsec) 


(K km s-^) (k 


m s"-*^) (ki 


m s"-*-) 


(K) 




CO 2-1 


230.537970 0.7 


33 


7.77 


7.59 


4.51 


1.62 


CSO 


C^^O 2-1 


219.560352 0.8 


35 


0.55 


7.16 


0.40 


1.28 


CSO 


HCO+ 3-2 267.557620 0.7 


29 


0.36 


7.18 


0.48 


0.70 


CSO 


N2H+ 3-2 


279.511701 0.7 


27 








< 0.13 


CSO 


N2D+ 3-2 


231.321635 0.7 


33 








< 0.15 


CSO 


N2H+ 1-0 


93.176258 0.5 


54 


0.59 


7.90 


0.23 


0.41 


FCRAO 


CS 2-1 


97.980953 0.5 


54 


0.42 


7.39 


0.83 


0.54 


FCRAO 



Table 6. The best fit model parameters of CB130-1 

Parameter One-dimensional Two-dimensional 

Stellar Temperature 3000 K 

Internal Luminosity 0.16 ± 0.005 L© 0.14 ± 0.005 Lq 
Inner radius 350 AU 

Inclination angle . . . 50-70° 

Disk scale height ... 0.05 ± 0.005 
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Table 7. Best-fitting Step Function Abundance Parameters 



Line 


^0 




td ipc) 


Id 


CO 


7.40 X 10" 


-5 


0.04 


30 




1.37 X 10- 


-7 


0.04 


30 


HCO+ 


7.00 X 10- 


-8 


0.04 


10 




1.00 X 10- 


-9 


0.04 


10 


DCO+ 


1.00 X 10" 


-9 


0.04 


10 


N2H+ 


5.00 X 10- 


10 


0.04 


0.1 


cs 


1.00 X 10" 


-9 


0.04 


3 


H2C0 


2.00 X 10" 


-8 


0.04 


10 



Note. — The free parameter values 
of the step function abundance profiles 
obtained with the line modeling. 
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Fig. 1.— The DSS R band image (gray scale) with CS (2-1) (left) and N2H+ (1-0) (right) 
contours from the FCRAO (De Vries et al., in preparation). The CS (2-1) contour levels are 
0.4, 0.6, 0.8 times the peak value of 0.32 K km s'^ The N2H+ (1-0) contour levels are 0.3, 
0.4, 0.5, 0.6, 0.7, 0.8, 0.9 times of the peak value of 1.37 K km s-\ N2H+ (1-0) emission is 
not detected from the CB 130-2 core. The beam size is shown in the lower right corner of 
both panels. 
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Fig. 2. — Three color image of CB130-1 using the cores2deeper IRAC image with IRAC 3.6 
fim, 4.5 /im, and 8.0 /im as blue, green and red, respectively. CB130-1-IRS1 and IRS2 are 
labeled and marked with black plus signs. 



ie''16"'19' 18"! e"! 8* IB^Ciy 18*1 6"1 6° 18*1 6"! 5" 

RA (J2000,0) 



Fig. 3. — The three color image composed of near IR band image, with J as blue, H as 
green, and K as red. The left plus sign indicates CB130-1-IRS2, and the right plus sign is 
the position of CB130-1-IRS1. The light at the position of CB130-1-IRS1 is mostly scattered 
light from CB130-1-IRS1. The white arrows indicate the cone shape nebulosity from CB130- 
1-IRSl extending toward CB130-1-IRS2. 
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Fig. 4. — The photometry images from the different wavelengths. The J, H, K band images 
are from ISPl, the 3.6 /xm, 4.5 /im, 5.8 /im, and 8.0 fim images are from IRAC, and 24 fim 
image is from MIPS. We plot the contours taken from the 850 fim SCUBA. The beam size 
of SCUBA at 850 /xm is shown in the bottom right of the 24 /xm panel. The contour levels 
are 0.5, 0.6, 0.7, 0.8, 0.9 times of the peak value (248 mJy/beam). 
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Fig. 5. — A color-color diagram for the 4 IRAC bands. The notation [3.6] — [4.5] denotes 
the magnitude difference, or color, between 3.6 and 4.5 fim. The box indicates the location 
of Class II sources (Allen et al. 2004). The two YSOs in CB130-1 are clearly redder than 
background stars. CB130-1-IRS2 is clearly in the box where Class II sources are found, 
whereas IRSl is not. The black dots are classified as stars (those with colors near zero) or 
galaxies (those in the Class II or Class I regions). 
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Fig. 6.— The SED of CB130-1-IRS1 (upper panel) and CB130-1-IRS2 (lower panel). The 
photometry data from ISPI, IRAC, MIPS, SHARC-II and SCUBA are plotted as filled circles 
with error bars. The IRS spectrum of CB130-1-IRS1 is plotted as a solid line. CB130-1-IRS1 
has a generally rising infrared SED, and peaks at 350 fim. CB130-1-1RS2 is not detected 
longward of 24 fira. The dashed line in the lower panel is the slope of the SED from 2.17 
/xm to 24yLim. The slope is —0.69, which indicates CB130-1-IRS2 is a Class II object. 
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Fig. 7. — The observed molecular lines of CB130-1 and line modeling results for step 
function abundances (see The observed data are plotted as solid lines. The vertical 
dashed lines are line centers determined by averaging the LSR velocity of optically thin 
lines. The modeled data are plotted as dotted lines. The modeled C^^O, H^^CO, N2H"'", and 
CS lines show a good fit with the observed lines, while CO and HCO"'" lines show a poor fit 
with the observed lines. CO is confused by the wing component; modeled lines of HCO"'" 
and H2CO show self- absorption features. 
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Fig. 8. — The observed molecular lines of CB130-2, and CB130-3. The parameters and 
descriptions are in Table |4] and Table |5| 
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Fig. 9. — The photometry data and the best fit 1-D radiative transfer model resuh for 
CB130-1-IRS1. The photometry data are plotted as filled circles. The data from the IRS is 
plotted as a solid line. The dash-dot line is the input SED from the internal source. The 
final SED is plotted as the dashed line. The aperture size convolved SED is marked with 
open circles at wavelengths longer than 100 /im. The model parameters are T^ = 3000 K, 
Lint — 0.16 Lq, Rin — 350 AU. The lower panel indicates the radiative transfer modehng 
result with a passive disk model. With a temperature distribution of a passive disk model, 
we cannot fit the 24 /im data. 
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Fig. 10. — The plots for a grid of one-dimensional radiative transfer models. Panel (a) is 
versus Tg, with L^t and Rint held fixed at 0.16 L© and 350 AU, respectively. The smallest 
is at Ts = 3000 K. Panel (b) is versus L^, T, and i^^n held fixed at 3000 K and 350 

AU. The smallest is at L^nt = 0.16 Lq. Panel (c) is versus Rin, with Tg and L^nt held 

fixed at 3000 K, and 0.16 Lq. The smallest is at R^n = 350 AU. 



- 39 - 



4O0O 
30001 

■^t 2O0O 
1O0O 



(°) H„Bi/R„^ = D.D5 








0,10 0,15 0,ZD 0,25 



60 . 
.SO -■ 

20 - 



v. M "H:LH.'.! "vlJ', 



30 40 50 60 70 
Inclination^degree) 



80 



Fig. 11. — The plots for a grid of two-dimensional radiative transfer models. Panel (a) is 
versus Hdisk/ Rdiski with Lint — 0.14 L0, and inclination angle 45 degrees. The inclination 

angle is face-on. The smallest is at Hdisk/Rdisk — 0.05. Panel (b) is versus Lint, with 
Hdisk/Rdisk = 0.05, and inclination angle 45 degrees. The smallest is at L^nt = 0.14 Lq. 
Panel (c) is versus inclination angle, with Hdisk/Rdisk = 0.05, and Li^t = 0.14 Lq. The 
smallest x^ is at inclination angle 70 degrees. 
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Fig. 12. — SEDs obtained from the 2-D radiative transfer modeling. The red filled circles 
are the photometry data from the Spitzer, CSO, and JCMT, and the red solid line is the 
spectrum from the Spitzer IRS. The black solid line is the SED obtained from the model at 
inclination angle of 50°. The blue open circles are flux densities at wavelengths where the 
source model has been convolved with the aperture used for measuring fluxes. The green 
dashed line shows the model SED for an inclination angle of 70° for comparison. 
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Fig. 13. — The first panel is the density profile as a function of radius, which follows 
n(r) oc r~^ ®. The second panel is the gas and dust temperature distribution as a function of 
radius. The solid line is the dust temperature and the dash-dot line is the gas temperature 
profile. The dust temperature and the gas temperature are almost the same in the inner 
region due to the high collision rate, because of the high density in the inner part. In the 
outer part of the cloud the gas temperature is lower than the dust temperature, since the 
photoelectric heating is too low to heat the gas. The last panel is the step function abundance 
profiles as a function of radius. The molecules freeze out in the inner part of the cloud. Since 
CO destroys the N2H"'" molecule, the abundance profile of N2H"'" has reversed shape to other 
molecules. 
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Fig. 14. — The luminosity evolution of the models we have used for chemical evolution. 



The dash-dot line is Model 1 (Young & Evans 2005), the dashed line is Model 2 (Lee et al. 



2004), and the solid line is Model 3 ( Dunham et al.|[2"010b ). The early stage evolution is the 



same in all three models. The two horizontal lines indicate the best fit luminosities for 1-D 
and 2-D models 
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Fig. 15. — The abundance profiles with Model 1, Model 2, Model 3, Model 4, and the step 
function abundances for comparison. The abundance profiles with Model 1 at 20000 years 
are plotted in solid lines. The abundance profiles from Model 2 at 36000 years are plotted in 
dash-dot lines. The abundance profiles from Model 3 at 78000 years are plotted in dashed 
lines. The dash-dot-dot-dot lines are the abundance profiles of Model 4. For comparison, 
we plotted the best fit step function abundances in thin dotted lines. The chemical modeled 
abundance profiles are in the range of the step function abundances except for HCO"^ and 
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Fig. 16. — The molecular lines from the chemical evolution model with Model 1. The solid 
lines are the observed molecular lines and the dotted lines are modeled molecular lines. The 
time is 20000 yr with luminosity 0.011 L0, which is just before the FHSC stage ends. 
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Fig. 17. — The molecular lines from Model 2. The solid lines are the observed molecular 
lines and the dotted lines are modeled molecular lines. This model includes a long FHSC 
transition stage. The time is 45000 yr with the internal luminosity 0.16 Lq. 
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Fig. 18. — The molecular line result at 78000 years with Model 3. The solid lines are 
observed molecular lines, and the dotted lines are the modeled lines at 78000 yr of Model 3. 
The internal luminosity at this stage is 0.14 L0. It has gone through six periods of accretion 
bursts, and 3000 years have passed after the accretion burst. The modeled C^^O is stronger 
and the modeled N2H"'" is weaker than the observed lines. 
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Fig. 19. — The molecular line result with Model 4. The solid lines are observed molecular 
lines, and the dotted lines are the modeled lines. The evolution stage is the same as Model 
3. The calculation turned 80% of CO ice into CO2 ice during low luminosity periods. The 
lines of C^^O, CS, DCO+, H2CO, and N2H+ are fitted better than in the other models. 



